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Abstract. 

A new, momentum preserving fast Poisson solver for N-body systems sharing ef- 
fective 0{N) computational complexity, has been recently developed by Dehnen 
(2000, 2002). We have implemented the proposed algorithms in a Fortran-90 code, 
and parallelized it by a domain decomposition using the MPI routines. The code 
has been applied to intensive numerical investigations of galaxy mergers, in par- 
ticular focusing on the possible origin of some of the observed scaling relations 
of elliptical galaxies. We found that the Fundamental Plane is preserved by an 
equal mass merging hierarchy, while it is not in a scenario where galaxies grow 
by accretion of smaller stellar systems. In addition, both the Faber-Jackson and 
Kormendy relations are not reproduced by our simulations. 
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1. Introduction 

In the last few years, fast algorithms for 
computing N-body interactions relying on 
the general multipole expansion techniques 
have been developed. These schemes, usu- 
ally referred to as Fast Multipole Methods 
(FMMs, see, e.g., Greengard & Rokhlin 
1987, 1997), are designed to have 0{N) 
computational complexity. However, it has 
been noted that only 0{NlogN) scaling is 
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usually achieved in numerical implementa- 
tion (see, e.g., Capuzzo-Dolcetta & Miocchi 
1998, hereafter CM98). The 0(iVlog7V) op- 
eration count also characterizes the N-body 
codes based on tree algorithms, originally 
proposed by Barnes & Hut (1986, hereafter 
BH86). In addition, it has been shown by 
direct numerical tests that, for given accu- 
racy, tree-codes are faster than the FMMs 
by a factor of 4 (CM98) . For this reason, 
most of the N-body codes currently used 
for astrophysical applications are based on 
the classical BH86 tree scheme, where the 
particle data are organized into a nested 
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oct-tree cell structure. Interactions between 
distant particles are approximated by cell- 
particle interactions, and the distribution 
of the particles on each cell is represented 
by a multipole expansion, usually trun- 
cated to the quadrupole term. 



Recently, a new scheme has been intro- 
duced (Dehncn 2000, 2002, hereafter D02), 
which can be seen as an original combina- 
tion of the tree-based BH86 scheme and 
of the FMMs multipole expansion, trun- 
cated at a fixed low order level (the pro- 
posed p = 3 order results the optimal 
one). This scheme, implemented by the 
Author in a C-|— |- code [named falcON, 
Force Algorithm with Complexity 0{N)], 
results in a significant improvement over 
the existing BH86 tree-codes, in terms of 
uniform accuracy, linear momentum con- 
servation, and CPU time performances. 
Moreover, the numerical tests presented in 
the D02 paper show that for an accuracy 
level e ~ 10~^, for the first time an effec- 
tive 0{N) scaling in operation count, as 
predicted by the FMM theory, is in fact 
obtained. For these reasons two of us (P.L. 
and C.N.) have implemented this scheme in 
a new Fortran-90 code (FVFPS, a Fortran 
Version of a Fast Poisson Solver), and 
have then parallelized it by using the MPI 
procedures. A completion of the code to 
handle gas-dynamics using Euler equations 
and the AMR (Adaptive Mesh Refinement) 
techniques is under development. Here we 
briefly describe the main characteristics 
of Dehnen's algorithm, and we present 
our FVFPS implementation and its paral- 
lelization. The achieved accuracy level, the 
performances in terms of operation count 
speed-up, the effective 0{N) scaling, as 
well as the parallel scaling with the proces- 
sors number, are then documented in nu- 
merical tests. Finally, we describe some in- 
teresting astrophysical applications of the 
code, in the context of the study of galaxy 
merging. Our version of the code is avail- 
able upon request. 



2. Dehnen's scheme 

In the FMM approach, a system of particles 
is first organized in a oct-tree structure of 
nested cells, providing a hierarchical group- 
ing of particles. For each pair of distant 
cells {A,B), the interactions between par- 
ticles are then approximated by a Taylor 
expansion of the two-point Green function 
around the geometrical center of each cell 
by using spherical harmonics. In these cell- 
cell interactions the two-body symmetry of 
the exact Green function is preserved by 
the multipole approximation. In principle, 
a FMM-based scheme can achieve a pre- 
scribed accuracy e by increasing the order 
p of the Taylor polynomial, and/or by low- 
ering the expansion parameter t] = r/R 
(where R is the distance between the two 
cell centers and r gives the typical size of 
a cell enclosing a particle group). Even if 
a formal, asymptotic scaling 0{N) is ex- 
pected, the proposed algorithms in the orig- 
inal formulation (Greengard & Rokhlin, 
1997) appear in fact to be less perform- 
ing and of difficult application. Therefore, a 
first variant has been introduced in CM98, 
where the expansion order is kept fixed to 
p = 2, and only the r/ parameter c;ontrols 
the accuracy of the Taylor expansion. This 
approach takes advantage of the adaptivity 
of tree algorithms to select the pairs of cells 
(even belonging to different levels of refine- 
ments), which satisfy the opening criterium 
rcTit/R < 0, where 9 is a preassigned con- 
trol parameter. To improve the efficiency of 
the scheme, a new definition of the critical 
size rcrit(j4, -B), has been also introduced. 
However, even in this improved formula- 
tion, the FMM scheme results to be not 
competitive with the classical BH86 tree- 
code. 

Dehnen's scheme basically develops 
along similar lines, but with some new in- 
gredients, which in fact succeeded to im- 
prove the overall eflaciency in a substantial 
way: 

— the Taylor polynomial expansion is per- 
formed in Cartesian components, and 
the expansion centers are the center of 
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mass of each cell, as in the BH86 tree- 
code; 

— a fixed expansion order p = 3 is chosen 
as optimal; 

— a new definition of the critical radius 
Tcrit is introduced for each cell, to opti- 
mize the adaptive selection of all pairs 
of well-separated cells; 

— a new, mass dependent control parame- 
ter = 9{M) has been introduced, that 
assures faster performances and a more 
uniform error distribution. The 6{M) 
function is given in implicit form by 



(1-^)2 (l-^„,in)2 VMtot; ' 

(1) 

where ^min (the minimum value of 9, as- 
sociated with the total mass) enters as 
a new, preassigned control parameter, 
and where M is the total mass enclosed 
in a cell. 

The resulting force solver is then organized 
as follows: 

(a) As in the BH86 code, an oct-tree 
structure of nested cells containing parti- 
cles is first constructed. A generic cell A, 
when acting as a gravity "source" is charac- 
terized by the position of its center of mass 
X^, by its critical size ta, mass Ma, and 
by the quadrupole tensor Qa- For the same 
cell A, when acting as a gravity "sink" a set 
of (twenty) Taylor coefficients Ca are used 
to store the potential and acceleration val- 
ues at X^. 

(b) To compute these coefficients, a new 
procedure for the tree exploration has been 
designed: starting from the root cell, all 
pairs of cells (A, B) satisfying the accep- 
tance condition {rA+rB)/R < 9, are looked 
for in a sequential way. When the condi- 
tion is satisfied, each cell of the selected 
pair accumulates the resulting interaction 
contributions in Ca{B) or Cb{A). In case 
of nearby cells (not satisfying the accep- 
tance condition), the bigger cell of the pair 
is splitted, and the new set of pairs are then 
considered. However, in order to avoid vis- 
iting the tree too deeply, some empirical 



conditions are adopted to truncate the se- 
quence, by computing directly (exact) two- 
body forces for particles in the smallest 
cells. The Author, by using a stack struc- 
ture to order the visited cell-pairs, has in- 
troduced a powerful and efficient new al- 
gorithm able to perform this interaction 
phase in the force computations. 

(c) Finally, once all the cell-cell interac- 
tions has been computed and stored in the 
Ca array, a final step is needed to evaluate 
the gravitational potential and the acceler- 
ations of the particles in the (generic) cell 
A. This evaluation phase can be performed 
by a simple and fast 0{N) tree-traversal. 

In our F-90 implementation of the 
scheme above, we followed rather closely 
the algorithms as described in D02: at this 
stage of development and testing, we have 
neither attempted to optimize our imple- 
mentation, nor to look for particular tricks 
to save computational time. 

To parallelize the force solver, we have 
adopted a straightforward strategy, by first 
decomposing the physical domain into a set 
of non-overlapping subdomains, each con- 
taining a similar number of particles (a few 
per cent of tolerance is allowed). The set 
of particles on each subdomain, Sk , is then 
assigned to a specific processing unit, Pk- 
Each processor builds its own tree [step 

(a) ], and computes the interaction phase 
for all the particles in its subdomain [step 

(b) ]. Steps (a) and (b) are executed in par- 
allel. By ordering the processors in a one- 
dimensional periodic chain, it is then possi- 
ble to exchange (particles and tree-nodes) 
data between them, by a "send_receive" 
MPI routine, in such a way that each pro- 
cessor can compute the interaction phase 
with the other subdomains. This phase is 
then performed by the following computa- 
tional steps: 

(i) The Pk processor sends its source- 
data to Pfe+s and receives from the Pk-s 
a copy of its source-data, where s = 
1, 2, ...A^pe/2, and A^pe is the total number 
of processors. 

[ii) At any given level s, each Pk com- 
putes, again in parallel, the interaction 
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phase with particles and cells of Pks- 
Mutuality of the interactions allows to eval- 
uate, at the same time, the results of these 
interactions for cells and particles of pro- 
cessor Pk-s to be stored as sink-data. 

{Hi) These accumulated sink-data (ac- 
celeration for particles and Taylor's coeffi- 
cients for cells) are sent back to the Pks 
processor, while the corresponding sink- 
data of the Pk processor are received from 
processor Pk+s and added to the current 
Pk values. In this way, it is evident that by 
repeating step {it) A^pe/2 times and step 
{Hi) Npe/2 — 1 times, the interactions of 
Pk particles with all other particles in the 
domain are recovered. 

{iv) The final evaluation phase [step 
(c)], being fully local, is performed in par- 
allel. 

Besides the force solver, the N-body 
code has been completed by a leap-frog 
time integrator scheme, with a (uniform) 
time step size At determined adaptivcly 
using a local stability condition At < 
1 / v'47rGpinax, where G is the gravitational 
constant, and pmax the maximum of par- 
ticle density. For typic;al simulations in- 
volving galaxy mergers, this adaptive (in 
time) step size keeps comparable to the 
value Tdyn/100, where T^yn represents a 
macroscopic dynamical time scale, and as- 
sures energy conservation with relative er- 
rors smaller than 0.1%. 

3. Performances of the FVFPS code 

We ran some test simulations in order to 
estimate the performances of our FVFPS 
code, and to compare it both with fal- 
cON and with the BH86 tree-code. We 
focus first on the analysis of the time 
efficiency in force calculation, consider- 
ing (one-processor) scalar simulations of 
the N-body system representing a galaxy 
group, where each galaxy is initialized by 
a Hernquist (1990, hereafter H90) galaxy 
model. Values of the total number of par- 
ticles A'' up to 2 X 10® are considered. The 
opening parameter given in eq. (1) is used 
by adopting 9min = 0.5, while for the BH86 



code we used 6 = 0.8: this choice results in 
comparable accuracies in the force evalua- 
tion. We ran our simulations on a Pentium 
III/1.133 GHz PC. For falcON and for the 
BH86 code we refer to the data published 
by D02, who used a Pentium III/933 MHz 
PC. 

In Fig. 1 (left panel) we plot the time 
efficiency of the codes as a function of the 
number of particles; it is apparent from the 
diagram that FVFPS (solid line), as well 
as falcON (dashed line), reaches an effec- 
tive 0{N) scaling, for Af>10l Our F-90 
implementation results to be slower than 
the C-I--I- version by a factor of ~ 1.5: we 
think that this discrepancy may be due to 
differences in the details of the implemen- 
tation, and also reflects a better intrinsic 
efficiency of the C-|— I- with respect to the F- 
90 for this kind of algorithms (especially for 
routines where intensive memory accessing 
is required). In the same diagram the data 
relative to the BH86 tree-code are also plot- 
ted (dotted line), and the 0{NlogN) scal- 
ing is apparent. As already documented in 
D02, it is remarkable how, for fixed A'' and 
prescribed accuracy, the BH86 tree-code is 
significantly slower than the new scheme: 
for example for N ~ 10^ the FVFPS code 
is faster by a factor of ^ 8. Therefore, even 
if the present F-90 implementation does 
no yet appear to be optimal, it gives fur- 
ther support to the efficiency of the D02 
scheme, independently of the adopted pro- 
gramming language and of the implemen- 
tation details. 

The accuracy in force evaluation is doc- 
umented on the right panel of Fig. 1, where 
we plot, as a function of TV, the mean 
relative error on the acceleration modulus 
(filled circles) 

/ \ = ^ ST^ kapprox ~ QdirectI /r,\ 

and the corresponding 99th percentile 
(Aa/a)99% (empty squares), for FVFPS 
(solid lines) and, for comparison, in fal- 
cON (dashed lines). As usual, the refer- 
ence acceleration values ctdirect evalu- 
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Fig. 1. Left panel: CPU time per particle as a function of total number of particles for 
BH86 tree-code (dotted line), falcON (dashed line; data from D02), and FVFPS (solid 
line). For the BH86 code 6 = 0.8, while in the last two cases 6 = 0{M), with ^min = 0.5. 
Right panel: Mean error on the force calculation (filled circles) and the 99th percentile 
(empty squares), as a function of the number of particles in falcON (dashed line; data 
from D02), and in FVFPS (solid line). 



ated by direct summation of the exact two- 
body Green function, while the approxi- 
mated values are given with 6 = 6{M), 
and ^min = 0.5. In this case, we have 
considered only simulations having a sin- 
gle H90 galaxy model. From the diagram 
it is apparent that, for both implementa- 
tions, (Aa/a) and (Aa/a)99% do not exceed 
a few 10^'^ and a fc;w pc;r cent, respectively, 
assuring accuracies usually considered suf- 
ficient for astrophysical applications (see, 
e.g., CM98). 

As outlined in the previous Section, we 
parallelized our code by using a domain de- 
composition with a load balancing essen- 
tially based only on particles number on 
each subdomain. This load balancing cri- 
terium seems in fact to work well, taking 
advantage of the intrinsic 0{N) scaling in 
the operation count of the basic scheme. 
In order to check the performance of the 
parallelized algorithm, we ran some simu- 
lations on the IBM Linux cluster (Pentium 
III/1.133 GHz PCs) at the CINECA cen- 
ter, in a range A^pe < 16 in the num- 
ber of processors. In this case the initial 



conditions are represented by a group of 
8 one-component H90 galaxy models with 
a total number of particles N = 1048576. 
As can be observed in Fig. 2 (left panel), 
this parallelization is characterized by a 
rather well behaved linear scaling of the 
CPU time per particle versus A^pe. This 
scaling can be better evaluated by consid- 
ering the right panel of Fig. 2, where dif- 
ferent contributions to the total CPU con- 
sumption are plotted. As expected, up to 
A'pE = 16 the main contributions c;ome 
by the tree building ("tree") and by the 
the force computation in each sub-domain 
("self"), corresponding to the step (i) and 
{ii), respectively, of the parallel scheme 
outlined in the previous Section. The con- 
tributions coming from the force evaluation 
among particles belonging to different sub- 
domains ("cross") and the cumulative time 
for all communication steps among proces- 
sors ("trans"), which increase with A^pe, 
keep smaller in the chosen range of parame- 
ters (A^, A^pe). Therefore, the 0{N) scaling 
of the base scheme, allows to estimate, for 




Fig. 2. Left panel: number of particle processed per second by the parallel version of the 
FVFPS code, as a function of the number of processors. The data refer to simulations of 
a group of galaxies with total number of particles N — 1048576. Right panel: CPU time 
per particle as a function of the number of processors for the same simulations in left 
panel. The different contributions are plotted (see text for an explanation). 



given TV, the optimal configuration of the 
iVpE that can be used. 

4. Simulations of galaxy merging 

As an astrophysical application of the de- 
scribed FVFPS code, we consider here the 
study of the galaxy merging process. We 
checked the reliability of the numerical sim- 
ulations performed with FVFPS for this 
kind of application (now time evolution is 
also considered), by running a few merg- 
ing events, between spherically symmet- 
ric stellar systems, both with GADGET 
(Springel, Yoshida & White 2001), and 
with FVFPS. The results are in remarkable 
agreement: the remnants obtained with the 
two codes are practically indistinguishable 
for three-dimensional shape, circularized 
half-mass radius, circularized density and 
velocity dispersion profiles. As an illustra- 
tive example, in Fig. 3 we show the density 
profiles of the end-product of the merger 
obtained from the head-on parabolic en- 
counter of two H90 models, with a total 
number of particles N = 65536 obtained 
with GADGET (empty symbols), and with 



FVFPS (filled symbols); note how also "mi- 
nor" details, as the slope change in the den- 
sity profile at r ~ lOrM, are nicely repro- 
duced. In addition, in both cases the change 
in total energy AE/E is smaller than 10"''' 
per dynamical time, over lOOTdyn. In this 
test, when using the FVFPS code, we 
adopted softening parameter e — 0.05, ini- 
tial time step At = 0.01, and (constant) 
9 = 0.65; in the GADGET simulation we 
used e — 0.05, atoi = 0.05, Atj^in = 0, 
Atmax — 0.03 , a — 0.02 (see Springel et 
al. 2001 for details), where the softening 
parameters and the time step parameters 
are in units of the galaxy core radius and 
dynamical time, respectively. 

On the basis of the results of these tests, 
we used our N-body code to investigate 
with numerical simulations the effects of hi- 
erarchies of dissipationless galaxy mergers 
on the scaling relations of elliptical galax- 
ies. This study is described in details in 
Nipoti, Londrillo & Ciotti (2002, hereafter 
NLC02). Here we report, as an example, 
the results relative to the effects of merg- 
ing on the Fundamental Plane of elliptical 
galaxies (hereafter FP; Djorgovski & Davis 



6 



P. Londrillo et al.: N-body simulations of galaxy mergers 



a 

















'XT 








"a 


-2 






o 





- ■ = FVFPS ■□ - 

_g 1 □ = GADGET . 1 

U I I M?nr 

0.01 0.1 1 10 100 

Fig. 3. Angle-average density profile of the 
remnant of the merging of two identical 
H90 galaxy models with FVFPS (full sym- 
bols), and with GADGET (empty sym- 
bols), tm is the half-mass radius. 

1987, Dressier et al. 1987). We considered 
a hierarchy of 5 equal mass mergers and 
a hierarchy of 19 accretion events, corre- 
sponding to an effective mass increase of 
a factor of ^ 24 and ~ 16, respectively: 
in the former case the successive gener- 
ations are produced by merging pairs of 
identical systems obtained by duplicating 
the end-product of the previous step, while 
in the latter case each end-product merges 
with a seed galaxy. Each merging is pro- 
duced in a head-on parabolic encounter. 
The seed galaxies are spherically symmetric 
one-component H90 galaxy models, with 
number of particles N = 16384 each. As a 
consequence of the merging hierarchy, the 
number of particles in the simulations in- 
creases with the galaxy mass: in the last 
merging of equal mass systems and of the 
accretion scenario, the total number of par- 
ticles involved is of the order of 5.2 x 10^ 
and 3.2 x lO'^, respectively. 

In Fig. 4 we plot the results of equal 
mass merging and accretion simulations 
in the plane (ki^k^), where the observed 
FP is seen almost edge-on, with best-fit 
fca = 0.15fci -|- 0.36 (solid line) and scatter 



* = equal mass merging 
X = accretion 




3 3.5 4 4.5 



Fig. 4. The merging end-products in the 
{ki,ky,) plane, where the solid line rep- 
resents the edge-on FP relation with its 
observed l-cr dispersion (dashed lines). 
Equal mass mergers are shown as trian- 
gles; crosses represent the accretion hier- 
archy. Bars show the amount of projection 
effects. 

rms(fc3) ~ 0.05 (dashed lines; see Bender, 
Burstein & Faber 1992). In the diagram the 
end-products of equal mass mergers and 
accretion are identified by triangles and 
crosses, respectively. The (common) pro- 
genitor of the merging hierarchies (point 
without bar) is placed on the edge-on FP. 
The straight lines in Fig. 4, associated with 
the merger remnants, represent the range 
spanned in the (/ci, k^) space by each end- 
product, when observed over the solid an- 
gle, as a consequence of projection effects. 
It is interesting to note that for all the end- 
products the projection effects are of the 
same order as the observed FP dispersion. 
In addition it is also apparent from Fig. 4 
how, in case of equal mass merging, kj, and 
fci increase with merging in a way full con- 
sistent with the observed FP tilt and thick- 
ness. On the other hand, in the case of ac- 
cretions, after few steps the end-products 
are characterized by a k^ decreasing for in- 
creasing fci, at variance with the FP slope 
and the trend shown by the end-products 



p. Londrillo et al.: N-body simulations of galaxy mergers 



7 



B 5 



1 1 1 1 1 
- i = equal mass 


1 1 1 1 1 1 

merging 


- X = accretion 


A - 




A 




A 


A 






A A 


A 






A 








X 

X 

X xX 








XX 








"^^^ 

x^ 




XXx^ 


M 




, , , , 1 





log M/M„ 

Fig. 5. Sersic best fit parameter m vs. total 
stellar mass of the end-products at stage i 
of the merging hierarchy. Symbols are the 
same as in Fig. 4 

of equal mass mergers. As a consequence, 
the last explored models (corresponding to 
an effective mass increase of a factor ~ 16) 
are found well outside the FP scatter. 

The different behavior of the end- 
products in the two scenarios is due mainly 
to effects of structural non-homology. This 
is clearly shown by an analysis of the 
projected stellar mass density profiles of 
the end-products. In fact, we fitted (over 
the radial range 0.1 < R/ {R)c Ss 4) the pro- 
jected density profiles of the end-products 
with the Sersic (1968) i?^/™ law, and in 
Fig. 5 we plot the best fit parameter m as 
a function of the total mass of the systems, 
for equal mass mergers (triangles) and ac- 
cretion (crosses). Clearly, the best-fitting 
TO depends on the relative orientation of 
the line-of-sight and of the end-products 
of the simulations: the two points for each 
value of the mass in Fig. 5 show the range 
of values spanned by to when projecting the 
final states along the shortest and longest 
axis of their inertia ellipsoids. In case of 
equal mass merging, the trend is to increas- 
ing with galaxy mass, as observed in real el- 
liptical galaxies (see, e.g., Berlin, Ciotti & 
Del Principe 2002 and references therein). 



On the contrary, in case of accretion the 
best Sersic to parameter decreases at in- 
creasing mass of the end-product, a behav- 
ior opposite to what is empirically found. 

As it is well known, elliptical galax- 
ies follow, besides the FP, also other scal- 
ing relations, such as the Faber-Jackson 
(hereafter FJ, Faber & Jackson 1976), and 
the Kormendy (1977) relations. It is re- 
markable that, in our simulations, both the 
equal mass merging and the accretion hier- 
archies fail at reproducing these two rela- 
tions. In fact, the end-products have cen- 
tral velocity dispersion lower, and effective 
radius larger than what predicted by the FJ 
and the Kormendy relations, respectively. 
In case of equal mass mergers, these two ef- 
fects curiously compensate, thus preserving 
the edge-on FP. In other words, the rem- 
nants of dissipationless merging like those 
described, are not similar to real elliptical 
galaxies, though in some cases they repro- 
duce the edge-on FP (see NLC02 for de- 
tails). 

5. Conclusions 

The main results of our work can be sum- 
marized as follows: 

— Our implementation of the D02 scheme 
in a F-90 N-body code (FVFPS) has ef- 
fective 0{N) complexity for number of 
particles iV>10'*, though its time effi- 
ciency is lower by a factor of ~ 1.5 than 
Dehnen's C++ implementation. 

— As shown by numerical tests, the paral- 
lelized version of the FVFPS code has 
good scaling with the number of proces- 
sors (for example for numbers of parti- 
cles N ~ 10^, at least up to 16 proces- 
sors). 

— We found very good agreement between 
the results of merging simulations per- 
formed with our FVFPS code and with 
the GADGET code. 

— From an astrophysical point of view, 
our high resolution simulations of hier- 
archies of dissipationless galaxy merg- 
ers indicate that equal mass merging 
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is compatible with the existence of the 
FP relation. On the other hand, when 
an accretion scenario is considered, the 
merger remnants deviate significantly 
from the observed edge-on FP (see also 
NLC02). 

— The different behavior with respect to 
the FP of the merger remnants in the 
two scenarios is a consequence of struc- 
tural non-homology, as shown by the 
analysis of their projected stellar den- 
sity profiles. 

— In any case, the end-products of our 
simulations fail to reproduce both the 
FJ and Kormcndy relations. In the case 
of equal mass merging, the combination 
of large effective radii and low central 
velocity dispersions maintains the rem- 
nants near the edge-on FP. 
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